#################################################
# Functions
#################################################

# Standardize variables
st_var <- function(x) {
  (x - mean(x, na.rm = TRUE)) / (2 * sd(x, na.rm = TRUE))
}

# Run event study models
run_ev_study <- function(data, damage_prefix, cluster="pro_com",weights=NULL) {
  
  variables <- paste0(damage_prefix, c("_04", "_09", "_19"))
  
  for (var in variables) {
    data[[var]] <- st_var(data[[var]])
  }
  
  formula_str <- paste0(paste("party_share ~", 
                              paste(variables, collapse = " + "), 
                              "+ as.factor(cod_prov):as.factor(anno)|pro_com+anno|0|+"),
                        cluster)
  formula <- as.formula(formula_str)
  
  model <- felm(formula, data = data, cmethod = 'cgm2', weights = weights)
  
  coefficients <- coef(model)
  conf_intervals <- as.data.frame(confint(model))
  
  measure  <- rep(damage_prefix, times = 4)
  year     <- c("2004", "2009", "2014", "2019")
  estimate <- c(coefficients[paste0(damage_prefix, "_04")],
                coefficients[paste0(damage_prefix, "_09")],
                0,
                coefficients[paste0(damage_prefix, "_19")])
  lower_ci <- c(conf_intervals[1, 1], conf_intervals[2, 1], NA, conf_intervals[3, 1])
  upper_ci <- c(conf_intervals[1, 2], conf_intervals[2, 2], NA, conf_intervals[3, 2])
  
  result <- data.frame(measure, year, estimate, lower_ci, upper_ci)
  
  return(result)
}

# Run models for three outcomes and create unified dataset
combined_results <- function(party) {
  
  result_combined <- NULL
  
  damage_prefixes <- c("damage_mean", "damage_max", "damage_binary")
  
  for (damage_prefix in damage_prefixes) {
    result <- run_ev_study(party, damage_prefix, "SLL_2011")
    
    if (is.null(result_combined)) {
      result_combined <- result
    } else {
      result_combined <- rbind(result_combined, result)
    }
  }
  
  return(result_combined) 
  
}

# Extract standard errors
se.coef.felm <- function(model){
  
  summary(model)$coefficient[,2]  
  
}

###################################
# R function for making a stacked coefficient plot from different regression estimates
# by Piero Stanig, based in part on Yu-Sung Su's coefplot code
# =====================================================================
# First version: November 22, 2006
# This version: Feb 2023

stackplot<- function(fit1, fit2=NULL, fit3=NULL, fit4=NULL,
                     fit5=NULL,fit6=NULL,
                     longnames=NULL, center.labels=FALSE,return.position=FALSE,
                     x.label="", y.label="", main.label="", 
                     sub.label="",
                     varnames=TRUE, estimators=FALSE,  
                     one.sd=FALSE,fifty=TRUE,
                     lwd.sd=1,max.n=15, offset=.1, smallish=.02,
                     force.scale=NULL,
                     box=FALSE, symmetric=TRUE, 
                     cex.var=1.2, cex.pts=1.5,  
                     suppress.intercept=TRUE, omit=NULL,
                     thin.box=FALSE,  cex.main=2,
                     restrict=NULL, suppress=NULL, ind.labels=NULL, 
                     pch.type=21, smallmargin=FALSE, compressed=FALSE, manyplots=FALSE
){
  
  if (class(fit1)[1]=="matrix"|class(fit1)[1]=="data.frame"){
    
    coeffs<-as.matrix(fit1)
    
    n.models<-ncol(coeffs)/2
    
    estimators<-FALSE
    
    
    if (is.null(longnames)){
      
      if (!is.null(row.names(fit1))){
        
        longnames<-row.names(fit1)
        
      }else longnames<-" "
      
    }
    
    n.vars<-nrow(coeffs)
    
    max.length<-n.vars
    
    if (max.length>40) {
      
      stop(message = "Can display up to 40 coefficients per model in a plot!\n")
      
    }
    
  }else{
    
    
    count.them<-c(is.null(fit1),is.null(fit2),is.null(fit3),is.null(fit4),is.null(fit5), is.null(fit6))
    
    n.models<-6-sum(as.numeric(count.them))
    
    if (n.models==1) fit<-list(fit1)
    
    if (n.models==2) fit<-list(fit1,fit2)
    
    if (n.models==3) fit<-list(fit1,fit2,fit3)
    
    if (n.models==4) fit<-list(fit1,fit2,fit3,fit4)
    
    if (n.models==5) fit<-list(fit1,fit2,fit3,fit4, fit5)
    
    if (n.models==6) fit<-list(fit1,fit2,fit3,fit4, fit5, fit6)
    
    for (k in 1:length(fit)){
      
      ### fix this with a "is.element"
      
      if (! (class(fit[[k]])[1]=="matrix"|class(fit[[k]])[1]=="data.frame"|
             class(fit[[k]])[1]=="felm")){
        
        stop(message = paste(class(fit[[k]]), "class not supported!\n"))
      }
    }
    
    n.predictors<-array(NA, c(1, n.models))
    
    coeffs<-array(NA, c(40,n.models*2))
    
    predictor.names<-array(NA, c(40, n.models))
    
    #### Extract the coefficients
    
    for (k in 1:n.models){
      
      object.class <- class(fit[[k]])[1]
      
      if(object.class=="logical"){
        
        coef.ml=rep(NA,3) 
        sd.ml =rep(NA,3)
      }
      
      if(object.class=="felm"){
        
        coef.ml<- summary(fit[[k]])$coefficients[,1]
        
        sd.ml <- summary(fit[[k]])$coefficients[,2]
        
      }else{
        
        coef.ml<- summary(fit[[k]])$coef[,1]
        
        sd.ml <- summary(fit[[k]])$coef[,2]
      }
      
      if(suppress.intercept==TRUE){
        
        coef.ml<- coef.ml[-1]
        
        sd.ml <- sd.ml[-1]
        
      }
      
      if(!is.null(omit)){
        
        coef.ml <- coef.ml[-omit]
        
        sd.ml <- sd.ml[-omit]
      }
      
      n.predictors[1,k]<-length(coef.ml)
      
      coeffs[1:length(coef.ml),(2*k)-1]<-coef.ml
      
      coeffs[1:length(coef.ml),2*k]<-sd.ml
      
      if(!is.null(names(coef.ml))){
        
        predictor.names[1:length(coef.ml),k]<-names(coef.ml)}
    }
    
    max.length<-max(n.predictors)
    
    n.vars<-max.length
    
    coeffs<-as.matrix(coeffs[1:max.length,])
    
    if(dim(coeffs)[2]==1){coeffs <- t(coeffs)}
    
  }
  
  if (max.length>40) {
    
    stop(message = "Can display up to 40 coefficients per model in a plot!\n")
    
  }
  
  ####Manage the variable names
  
  if (is.null(longnames)){
    
    longnames<-predictor.names[1:max.length,]
    
    longnames<-t(na.exclude(t(longnames)))[,1]
    
  }
  
  ####Stack the coefficients and the sd's
  
  coef.mat<-array(NA, c(n.models*n.vars,2))
  
  suppress.mat<-array(NA, c(n.models*n.vars,1))
  
  for (i in 1:n.models){
    
    for (j in 1:n.vars){
      
      coef.mat[((j-1)*n.models)+i,  1] <-coeffs[j,(i*2)-1]
      
      suppress.mat[((j-1)*n.models)+i,  1] <-is.element(j,suppress)
      
      coef.mat [i+(j-1)*n.models, 2]<-    coeffs[j,i*2]
      
    }
  }
  
  coef<-coef.mat[,1]
  
  sd<-coef.mat[,2]
  
  omit.from.plot <- as.numeric(suppress.mat)
  
  missing.c<-is.na(coef)
  
  missing.s<-is.na(sd)
  
  coef[missing.c==TRUE]<-0
  
  sd[missing.s==TRUE]<-0
  
  ###now stack them
  
  n.x <- length(coef)
  
  ####GROUP THE idx SO THAT THE COEFFICIENTS 
  ####ARE CLOSE TO EACH OTHER
  
  idx <- rep(seq(1 , n.x/n.models),each=n.models)
  
  smallsteps<-rep(.4*seq(0,1,length.out=n.models), n.x/n.models) 
  
  if(compressed==TRUE){
    
    idx<-idx+smallsteps*.5}else{
      
      idx<-idx+smallsteps}
  
  min.x <- round(min(coef-2*sd), 1)-0.1           
  
  max.x <- round(max(coef+2*sd), 1)+0.1                                 
  
  x.scale <- c(min.x-offset, max.x+offset)                                            
  
  y.scale <- range(idx)+c(0,.02)                   
  
  if(offset!=.1){symmetric=FALSE}
  
  if (symmetric) x.scale<-c(-max(abs(min.x), abs(max.x)), max(abs(min.x), abs(max.x)) )
  
  if(!is.null(force.scale)){ x.scale <- force.scale}
  
  subtitle<-""
  
  if (subtitle=="") subtitle<-sub.label
  
  notshown <- as.numeric(missing.c)+as.numeric(omit.from.plot)
  
  notshown[notshown==2]<-1
  
  if(manyplots==FALSE){
    
    if(smallmargin==FALSE){
      
      par(mar=c(4+1*as.numeric(estimators==FALSE), 
                2, 2+1*as.numeric(main.label!=""), 
                .2+.6*max(nchar(longnames)))   )
      
    }else{ 
      
      par(mar=c(3.5+1*as.numeric(estimators==FALSE), 
                
                2, 2+1*as.numeric(main.label!=""), 
                .1+.22*max(nchar(longnames)))   )  }
    
  }        
  
  plot(coef, idx, axes=FALSE, type="n",                                     
       
       xlim=x.scale, ylim=y.scale+c(-.02, .02), 
       
       xlab=x.label, ylab=y.label,                   
       
       sub=subtitle)                                                        
  
  title(main=main.label,cex.main=cex.main)
  
  axis(1, cex.axis=1.3)    
  
  abline(v=0, lty=2, col="grey40",lwd=1-.4*thin.box)                                                 
  
  if(n.models==1|length(pch.type)==1){
    
    points(coef, idx, pch=pch.type, cex=cex.pts*(1-notshown))
    
  } else{
    
    temp.i =1 
    
    for( i.dots  in 1:length(idx)){
      
      points(coef[i.dots], idx[i.dots], pch=pch.type[temp.i], 
             cex=cex.pts*(1-notshown[i.dots]))
      
      temp.i=temp.i+1 
      
      if(temp.i> length(pch.type)){temp.i=1}
      
    }
  }
  
  if (box==TRUE){
    
    for (k in 1:n.vars){
      
      abline(h=k-.2, lty=3, lwd=1-.4*thin.box)
    }
    
  }
  
  
  
  if(fifty==TRUE){crit=1}else{crit=1.645}
  
  if (one.sd==FALSE) segments ((1-notshown)*(coef+2*sd), 
                               idx, (1-notshown)*(coef-2*sd), idx, col="Grey80", lwd=lwd.sd)
  
  
  segments ((1-notshown)*(coef+crit*sd), idx, 
            (1-notshown)*(coef-crit*sd), idx, lwd=lwd.sd)     
  
  if(center.labels==TRUE){adj.here=.5}else{adj.here=1}
  
  if (varnames==TRUE) {
    
    axis(4, .2*as.numeric(n.models>1)+(1:n.x), 
         longnames[1:n.x], las=2, tck=TRUE, 
         pos=1.1*max(coef+2*sd), lty=0, adj=adj.here, 
         cex.axis=cex.var)
    
    out.now <- 1:n.x
  }
  
  if(!is.null(ind.labels)){
    
    text(coef, idx-smallish, ind.labels, cex=.8)  
    
  }		
  
  par(mar=c(5, 4, 4, 2) + 0.1)
  
  return(idx)
}